Theoretical principles underlying 1 optical stimulation of a 2 channelrhodopsin - 2 positive 3 pyramidal neuron 4 5 6 7 8 9 10 11

34 Optogenetics is an emerging field of neuromodulation that permits scaled, 35 millisecond temporal control of the membrane dynamics of genetically targeted 36 cells using light. Optogenetic technology has revolutionized neuroscience 37 research; however, numerous biophysical questions remain on the optical and 38 neuronal factors impacting the modulation of neural activity with photon sensitive 39 ion channels. To begin to address such questions, we developed a 40 computational tool to explore the underlying principles of optogenetic neural 41 stimulation. This “light-neuron” model consists of theoretical representations of 42 the light dynamics generated by a fiber optic in brain tissue, coupled to a multi43 compartment cable model of a cortical pyramidal neuron embedded with 44 channelrhodopsin-2 (ChR2) membrane dynamics. Simulations revealed that the 45 large energies required to generate an action potential are primarily due to the 46 limited conductivity of ChR2, and that the major determinants of stimulation 47 threshold are the surface area of illuminated cell membrane and proximity to the 48 light source. Our results predict that the activation threshold is sensitive to many 49 of the properties of ChR2 (density, conductivity and kinetics), tissue medium 50 (scattering and absorbance), and the fiber optic light source (diameter and 51 numerical aperture). We also illustrate the impact of redistributing the ChR2 52 expression density (uniform versus non-uniform) on the activation threshold. The 53 model system developed in this study represents a scientific instrument to 54 characterize the effects of optogenetic neuromodulation, as well as an 55 engineering design tool to help guide future development of optogenetic 56 technology. 57 Introduction 58 Optical stimulation technology has rapidly advanced since the first 59 characterization of Channelrhodopsin-2 (ChR2) (Nagel et al., 2003). ChR2 is a light 60 activated, nonspecific cation channel (H, Na, K, Ca) (Ehlenbeck et al., 2002; Zhang 61 et al., 2007; Berndt et al., 2010) that is now being expressed in vivo in a range of 62 mammalian species and can be targeted to specific neuron types. ChR2 photoactivation 63 begins with the absorption of blue spectrum light, followed by excitation of the retinal 64 chromophore, leading to a conformational change with opening of the ion channel 65 (Hegemann et al., 2005, Muller et al. 2011). Selective expression of ChR2 channels in 66 neurons and their activation by targeted delivery of light, permits millisecond temporal 67 activation of specific neural populations. This technique has been useful in elucidating 68 numerous sensory pathways, and providing new insights into the function of central 69 nervous system circuits (Deisseroth et al., 2006). 70 The scientific value of optogenetic technology is unquestioned; however, 71 numerous biophysical questions remain on the mechanism of action potential 72 generation from optical stimulation. Further, definitive relationships that describe light 73 irradiance thresholds as a function of the stimulation parameters are lacking on topics 74 as fundamental as the fiber-optic-to-neuron distance, the ChR2 75 density/conductivity/distribution, and the optical fiber geometry. In this study, we 76 investigate the complex relationship between an implanted optical fiber and a pyramidal 77 neuron using a “light-neuron” model. Our theoretical analysis attempts to expand 78 current knowledge, parameterizing the model based on the available experimental data, 79 and provides new quantitative hypotheses on the spatial extent of neural activation 80 induced by optogenetic stimulation. 81 Characterizing the underlying mechanisms of optical stimulation also has 82 important implications for the creation of new optogenetic technology. Knowledge of 83 neuron orientation and/or photon sensitive ion channel characteristics that best interact 84 with the light source to minimize stimulation power have important implications for the 85 engineering of next generation optical fibers and/or optrodes (e.g. Sparta et al., 2011), 86 as well as new channel constructs (e.g. Berndt et al., 2011). Further down the 87 optogenetic technology development line, proposals of therapeutic application in human 88 conditions such as Parkinson’s disease, spinal cord injury, depression and obsessive89 compulsive disorder will need to compete with the current neuromodulation 90 technologies which focus on electrical stimulation (Henderson et al., 2009). Numerous 91 advances in electrical stimulation technology can be attributed to the mechanistic 92 understanding provided by “field-neuron” models (McNeal, 1976). We propose that 93 “light-neuron” models could provide a similar service to the optogenetic community. 94 Preliminary portions of this work have been presented in abstract form (Foutz and 95 McIntyre, 2010; 2011). 96 Methods 97 Neuron model 98 The neuron model used in this study was based on the soma-dendritic cable 99 model of a reconstructed neocortical, layer V pyramidal neuron from cat visual cortex 100 (Mainen et al., 1995; Mainen and Sejnowski, 1996; Shu et al., 2006; Yu et al., 2008; Hu 101 et al., 2009) (Figure 1bc). The ModelDB accession number for the reconstructed 102 pyramidal neuron and endogenous membrane dynamics used in this study is 123897. 103 The neuron has an elliptical cell body (20 μm diameter), with 1 main apical dendrite and 104 tuft (1100μm in length), as well as a basal dendritic tree (Figure 1b). The main axon is 105 composed of an axon hillock (length, 10 μm; diameter, 3.8 to 2.4 μm), axon initial 106 segment (length, 50 μm; diameter, 1.22 μm), unmyelinated axon (length, 400 μm; 107 diameter 1.01μm) and myelinated axon (length, 1300μm; diameter, 1.21 μm; internodal 108 length, 100 μm). The myelinated axon incorporates 14 node compartments (length, 1 109 μm; diameter 0.91 μm). The membrane electrical properties were uniformly distributed 110 through the cell with the same values as previously published versions of the model (Hu 111 et al., 2009): specific membrane resistivity (rm) of 30 kΩ·cm, specific cytoplasmic (axial) 112 resistivity (ri) of 150 Ω·cm, and membrane capacitance (Cm) of either 1 μF·cm (soma), 113 0.5 μF·cm (dendrites, axon hillock, axon initial segment and nodes) or 0.02 μF·cm 114 (myelinated sections). The input resistance at the soma was 43.6 MΩ. Simulations were 115 run with a nominal temperature of 37 °C. 116 The endogenous membrane properties used in this study were identical to the 117 model reported in Hu et al. (2009), with the addition of exogenous ChR2 (Figure 1e; see 118 below) distributed uniformly throughout the neural structure (unless otherwise noted in 119 the results and figures). The neuron’s resting potential was -70 mV. Stimulation 120 threshold was defined using the irradiance exiting the tip of the fiber optic (see below). 121 The minimum irradiance necessary to generate an action potential that propagated to 122 the penultimate node of the axon was determined by a binary search algorithm. 123 Simulations were performed with NEURON v7.2 in Python (Carnevale and Hines, 124 2009). 125 126 Channelrhodopsin-2 dynamics 127 Channelrhodopsin-2 was modeled as a nonspecific ion channel with four states: 128 two closed states (C1, C2), and two open, conducting states (O1, O2) (Nagel et al., 129 2003; Hegemann et al., 2005; Nikolic et al., 2009) (Figure 1d). In this model, ChR2 can 130 be excited from a closed, ground state (C1) to an open state (O1) secondary to 131 absorption of a photon of ~470 nm light (Figure 1a). This process occurs with a rate 132 constant Ka1. ChR2 in the excited state (O1) can decay back to a closed state (C1, rate 133 constant Kd1) or transition to a second excited state (O2, rate constant e12). ChR2 in this 134 second excited state is more stable, but has a lower ion conductance. ChR2 in state O2 135 can either return to the first open state (O1, rate constant e21), or decay to a closed 136 state (C2, rate constant Kd2). Finally, channels in state C2 can either be photoexcited 137 back to O2 (rate constant Ka2), or be slowly converted thermally to C1 (rate constant Kr) 138 (Grossman et al., 2011). The instantaneous rate of change of these states was defined 139 by a set of rate equations: 140 dO1 dt = Ka1C1− (Kd1 + e12 )O1+ e21O2 (1) 141 dO2 dt = Ka2C2 + e12O1− (Kd2 + e21)O2 (2) 142 dC2 dt = Kd2O2 − (Ka2 + Kr )C2 (3) 143 1= O1+O2 +C1+C2 (4) 144 In these equations, O1, O2, C1 and C2 represent the fraction of ChR2 molecules in the 145 respective states. The fixed rate constants are summarized in Table 1. The activation 146 rate constants Ka1 and Ka2 are calculated dynamically, since they depend upon the light 147 irradiance. 148 Ka1 = ε1Φ(1− e − t /τ ), Φ > 0 ε1Φ0 (e −(t−t0 )/τ − e t /τ ), Φ = 0     (5) 149 Ka2 = ε2Φ(1− e −t /τ ), Φ > 0 ε2Φ0 (e −(t−t0 )/τ − e /τ ), Φ = 0     (6) 150 In these equations, ε1 and ε2 are the quantum efficiency of photons which attempt to 151 excite channelrhodopsin from a closed state to the corresponding open state; Φ is the 152 photon flux per unit area during illumination; Φ0 is equal to Φ during prior illumination 153 (dark phase); t is the time since prior illumination began; t0 is the time since prior 154 illumination ended (dark phase); τ is the time constant of channelrhodopsin, shown in 155 Table 1. 156 Determination of the ChR2 transmembrane channel conductance depends on 157 the transmembrane voltage (Vm), the reversal potential (Ecat, set to 0 mV), and the 158 channel conductance (gChR2). The ChR2 current during illumination (imax) is determined 159 by Ohm’s law: 160 imax = (Vm − Ecat )gChR2 (7) 161 ChR2 channel conductance is dependent on the state of the channel, with zero 162 conductance in states C1 and C2, low conductance (g2) in state O2, and high 163 conductance (g1) in state O1 (See Table 1). 164 After a period of illumination, the ChR2 transmembrane current decays 165 exponentially. This decay has been fit experimentally by separating the current into a 166 fast (ifast) and a slow component (islow) (Nikolic et al., 2009; Grossman et al., 2011). The 167 ChR2 transmembrane current, post-illumination, is defined as: 168 i = islowe −Λ1(t−t0 ) + i faste −Λ2 (t−t0 ) (8) 169 where Λ1 and Λ2 are current decay factors. As time increases, the net transmembrane 170 current decays to zero. The fast and slow components of the current are defined by: 171 i fast = imax O10 (Kd1 + (1−γ )e12,dark − Λ1)+Ο20 (γ (Kd2 − Λ1)− (1−γ )e21,dark ) Λ2 − Λ1 (9) 172 islow = imax Ο10 (Λ2 − Kd1 − (1−γ )e12,dark )+Ο20 ((1−γ )e21,dark +γ (Λ2 − Kd2 )) Λ2 − Λ1 (10) 173 where O10 and O20 are the fraction of open channels during the prior illumination 174 phase, and γ is the ratio of the conductance of the two states O2 and O1 (γ = g2/g1). 175 The ChR2 model was typically inserted in all compartments of the neuron model, 176 with a uniform ChR2 channel density of 130 μm (Nagel et al., 1995). To simulate the 177 effect of non-uniform distributions, ChR2 was distributed either in specific compartments 178 (Table 2), or by distributing its density based on distance from the soma (Figure 10b). 179 Distance-based ChR2 distribution was performed by weighting the channel density by 180 the path distance from the center of the soma to each point on the dendritic arbor, and 181 scaling the density such that the total number of the channels in the soma-dendritic 182 arbor remained constant. For the apical distribution, the most distant compartment had 183 the maximal density, while the soma had minimal density of ChR2. For the basal 184 distribution, the distribution was reversed. 185 186 Light Model 187 Most of our simulations were performed with the optical fiber oriented 188 perpendicular to the long axis of the neuron (Figure 1b), directed at the soma from a 189 distance of 1 mm, unless otherwise noted. There are four primary factors affecting the 190 distribution of light exiting the fiber optic. These are 1) the source light distribution, 2) 191 the geometric spread of unfocused light, 3) the scattering and 4) the absorbance of light 192 by the tissue (Figure 2). The light at each point in space (I) is defined by the source light 193 irradiance (I0, center of fiber optic output) and the transmittance of light between that 194 point and the source (T): 195 I(r, z) = T (r, z)I0 (11) 196 where r is the radial distance and z is the height in a cylindrical coordinate system with 197 the origin defined at the center of the fiber optic output. The transmittance is wavelength 198 dependent, and can be broken down into linear components: 199 T (r,z) =G(r,z)C(z)M (r,z) (12) 200 where G describes the Gaussian distribution of light emitted by fiber optics, C describes 201 the conical spreading of unfocused light, and M describes the scattering and 202 absorbance of light according to the Kubelka-Munk general theory of light propagation 203 (Kubelka and Munk, 1931; Vo-Dinh, 2003; Aravanis et al., 2007). The light-model 204 parameters are summarized in Table 1. 205 Light emitted from a fiber optic spreads as a cone of light (Figure 2a) with a 206 divergence half-angle (θdiv) dependent on the tissue index of refraction (ntis) and the 207 numerical aperture of the fiber optic (NAfo): 208 θdiv = sin −1 NAfo ntis     (13) 209 The radius of the light cone (R) at height z emitted by a fiber optic with radius R0 210 spreads according to: 211 R(z) = R0 + z tan(θdiv ) (14) 212 As the light diverges, the irradiance decreases according to the law of conservation of 213 energy. Therefore, when considering the effects of geometry independently, the radiant 214 power (P) is constant at all distances, and is equal to the irradiance (I) times the surface 215 area illuminated: 216 P = I(z)πR(z) = I0πR0 2 (15) 217 where I is the irradiance at distance z from the fiber optic. Therefore, the transmittance 218 due to geometrical spreading (C) is: 219 C(z) = R0 R(z)     2 (16) 220 The Gaussian distribution of light (G) emitted by a fiber optic (Weik, 1997) can be 221 approximated as a transmittance: 222 G(r, z) = 1 2π exp −2 r R(z)     2       (17) 223 According to this equation, 95.4% (2σ) of light is emitted by the core of the fiber optic, 224 with the remaining 4.6% emitted by the cladding. The Gaussian light distribution with 225 and without geometrical spread is shown in Figure 2a. 226 The last two factors affecting the light distribution in our model are the scattering 227 and absorptive properties of tissue. To capture these effects, we implemented the 228 Kubelka-Munk general model of light propagation in diffuse scattering media (Vo-Dinh, 229 2003). The transmittance of light in absorptive, scattering media (M) was: 230 M (r, z) = b asinh(bS r + z )+ bcosh(bS r + z ) (18) 231